function [results_fnames source_recon_results]=osl_run_source_recon_sensorspace(oat)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setup sensor space data
%
% results_fnames=osl_run_source_recon_sensorspace(oat)
%
% Mark Woolrich 2012

global OSLDIR;

dirname=oat.source_recon.dirname;
modalities=oat.source_recon.modalities;  % added by DM

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set first level diagnostic report up    
report_dir=[oat.results.plotsdir '/' oat.results.date '_source_recon'];
source_recon_report=osl_report_setup(report_dir,['Source recon (epoched) - sensor space data setup']);   

for sessi_todo=1:length(oat.source_recon.sessions_to_do),   

    sessi=oat.source_recon.sessions_to_do(sessi_todo);

    disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
    disp(['%%%%%%%%%%%%%%%%%%%%%%%  RUNNING OAT SOURCE RECON (SENSOR SPACE SETUP) ON SESS = ' num2str(sessi) '  %%%%%%%%%%%%%%%%%%%%%%%'])
    
    % set session specific diagnostic report up   
    report_dir=[source_recon_report.dir '/sess' num2str(sessi)];
    report=osl_report_setup(report_dir,['Session ' num2str(sessi)]);       

    source_recon=oat.source_recon;

    source_recon_sess=source_recon;
    
    source_recon_sess.do_plots=oat.do_plots;
    source_recon_sess.session_name=['session' num2str(sessi)];

    if ~isempty(source_recon.D_continuous),        
        [p fname e] = fileparts(source_recon.D_continuous{sessi});       
        source_recon_sess.D_continuous=[p '/' fname '.mat'];       
        disp('Using continuous data as input');        
    else
        source_recon_sess.D_continuous=[];
    end;
    
    if ~isempty(source_recon.D_epoched),
        [p fname e] = fileparts(source_recon.D_epoched{sessi});       
        source_recon_sess.D_epoched=[p '/' fname '.mat'];       
        disp('Using epoched data as input');
    else
        source_recon_sess.D_epoched=[];
    end;    
    source_recon_sess.mri=source_recon.mri{sessi};
    
    if length(source_recon.pca_dim)>1,
        source_recon_sess.pca_dim=source_recon.pca_dim(sessi);
    else
        source_recon_sess.pca_dim=source_recon.pca_dim;        
    end;
    
    clear source_recon;
            
    source_recon_results.source_recon=source_recon_sess;         
    source_recon_results.recon_method=source_recon_sess.method;

    % check SPM MEEG objects passed in
    only_epoched_data_provided=0;
    if ~isempty(source_recon_sess.D_epoched) && isempty(source_recon_sess.D_continuous),
        only_epoched_data_provided=1;
        source_recon_sess.D=source_recon_sess.D_epoched;
    elseif ~isempty(source_recon_sess.D_continuous),
        source_recon_sess.D=source_recon_sess.D_continuous;
    else
        error('Need to specify D_continuous or D_epoched.');
    end;

    % do epoching if needed, but only if any epoch info is provided
    do_epoching=0;
    if ~only_epoched_data_provided
        if ~isempty(source_recon_sess.D_epoched) || ~isempty(source_recon_sess.epochinfo),
            do_epoching=1;
        end;
    end;
    
    disp(['Running sensor space analysis for ' source_recon_sess.D]);
    disp(['Will be designated ' source_recon_sess.session_name]);
        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% move copy of SPM MEEG object to work on into the working directory
    warning off;
    mkdir(dirname);
    warning on;
    
    S2=[];

    S2.D=source_recon_sess.D;
    S2.newname=[source_recon_sess.session_name '_spm_meeg'];
    S2.updatehistory=0;
    D = spm_eeg_copy(S2);
    runcmd(['mv ' D.path '/' D.fname ' ' dirname]);
    runcmd(['mv ' D.path '/' D.fnamedat ' ' dirname]);
    D = spm_eeg_copy(S2);
    spm_filename=[dirname '/' D.fname];
    D=spm_eeg_load(spm_filename);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% do bandpass filtering
    S=source_recon_sess;
    if(length(S.freq_range)>1 && S.freq_range(2)>S.freq_range(1)),
        disp('Temporal filtering...');
        S2=[];
        S2.D=spm_filename;
        S2.use_fft_bandpass=1;

        if(0)

            % implement bandpass as a low pass followed by a high pass
            S2.filter.PHz=S.freq_range(2);
            S2.filter.band='low';
            Dnew = spm_eeg_filter_v2(S2);

            D.delete;
            D=Dnew;

            S2.D=D;
            S2.filter.PHz=S.freq_range(1);
            S2.filter.band='high';
            Dnew = spm_eeg_filter_v2(S2);

        else

            S2.filter.PHz=S.freq_range;
            S2.filter.band='bandpass';
            S2.use_fft_bandpass=1;
            Dnew = spm_eeg_filter_v2(S2);

        end;

        if S.bandstop_filter_mains
            % do notch filtering to remove mains noise
            notchBands = [48 52; 98 102; 148 152; 198 202; 248 252; 298 302; 348 352; 398 402; 448 452; 498 502];
            lowIn      = notchBands(:,1) > S.freq_range(1);
            highIn     = notchBands(:,2) < S.freq_range(2);
            doNotch    = find(lowIn & highIn);

            for iNotch = 1:numel(doNotch)
                S3              = [];
                S3.D            = Dnew;
                S3.filter.PHz   = notchBands(doNotch(iNotch),:);
                S3.filter.dir   = 'twopass';
                S3.filter.band  = 'stop';
                disp(['Band-stop filtering from ' num2str(S3.filter.PHz(1)) 'Hz to ' num2str(S3.filter.PHz(2)) 'Hz.']);
                Dnew = spm_eeg_filter_v2(S3);
            end % for iNotch = 1:numel(doNotch)

        end % if S.notch_filter_mains


        D.delete;
        D=Dnew;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% do epoching
    % epochinfo will be in source_recon_sess.epochinfo, if its empty then look in D_epoched for D_epoched.epochinfo

    if(do_epoching),
        S2 = [];
        S2.D = D;

        if(~isempty(source_recon_sess.epochinfo)),
            S2.epochinfo=source_recon_sess.epochinfo; 
        else
            try,
                D_epoched=spm_eeg_load(source_recon_sess.D_epoched);
                S2.epochinfo=D_epoched.epochinfo;
            catch
                error('No epoch info available');
            end;
        end;
        S2.epochinfo.padding = 0;
        S2.save=0;
        S2.reviewtrials=0;
        S2.bc=0; 
        disp('Doing no within-trial baseline correction at the point of epoching');

        Dnew = spm_eeg_epochs(S2);

        D.delete;
        D=Dnew;

        %% get bad channels and trials from passed in source_recon_sess.D_epoched    
        if(~isempty(source_recon_sess.D_epoched)),
            D_epoched_passed_in=spm_eeg_load(source_recon_sess.D_epoched);
            D = reject(D, 1:length(D.conditions), D_epoched.reject);              
            if(length(D_epoched_passed_in.badchannels)>0),
                D = badchannels(D, D_epoched_passed_in.badchannels, ones(length(D_epoched_passed_in.badchannels),1));
            end;
            D.save;
        end;
    end;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Select time windows of interest
    if isempty(source_recon_sess.time_range), 
        source_recon_sess.time_range = [D.time(1) D.time(end)];
    else
        if source_recon_sess.time_range(1)<D.time(1),
            error('source_recon_sess.time_range(1)<D.time(1)');
        end;
        if source_recon_sess.time_range(2)>D.time(end),
            error('source_recon_sess.time_range(2)<D.time(end)');
        end;
    end

    Events = D.events; GoodEpochs=[D.time(1) D.time(end)];
    if D.ntrials==1 % if continuous data then look for bad epochs
        for ev = 1:numel(Events)
            if isfield(Events,'type') && strcmp(Events(ev).type,'BadEpoch')
                GoodEpochs(end,2)=Events(ev).time;
                GoodEpochs(end+1,1)=Events(ev).time+Events(ev).duration;
                GoodEpochs(end,2)=D.time(end);
            end
        end
    end

    good_samples=zeros(1,D.nsamples);
    for i=1:size(GoodEpochs,1)
        good_samples(D.indsample(GoodEpochs(i, 1)):D.indsample(GoodEpochs(i, 2)))=1;
    end

    samples_of_interest=zeros(1,D.nsamples);
    for i=1:size(source_recon_sess.time_range,1)
        samples_of_interest(D.indsample(source_recon_sess.time_range(i, 1)):D.indsample(source_recon_sess.time_range(i, 2)))=1;
    end

    samples2use = samples_of_interest & good_samples;
    woi=[D.time(find(diff([0 samples2use])==1))' D.time(find(diff([samples2use 0])==-1))'];

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Establish trials
    trials = 1:D.ntrials;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Normalise modalities using smallest eigenvalues 
    %% calculated using good channels and good trials, and over all woi  
    if not(strcmp(modalities{1},'EEG'))   % added by DM

        disp('Normalising modalities...');

        S2=source_recon_sess;
        S2.D = D;
        S2.samples2use=samples2use;
        S2.trials=trials;       
        S2.do_plots=1;
        
        [ Dnew pcadim tmp norm_vec normalisation fig_handles fig_names] = normalise_sensor_data( S2 );
          
        % diagnostic plot of design matrix    
        report=osl_report_set_figs(report,fig_names,fig_handles);
        report=osl_report_print_figs(report);
    
        D.delete;
        D=Dnew;

    else
        normalisation=nan;
    end;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% HMM (do not use it in sensor space)

    NK=1;
    hmm_class=zeros(NK,size(D,2),size(D,3));
    hmm_class(1,find(samples2use),trials)=1;
    hmm_class_probs=zeros(NK,size(D,2),size(D,3));    
    hmm_class_probs(:,find(samples2use),trials)=1;
    
    % add channel
    Sc=[];
    Sc.D=D;
    Sc.newchandata=[hmm_class; hmm_class_probs];
    Sc.newchanlabels{1}='Class';
    Sc.newchantype{1}='CLASS';

    for kk=1:NK,
       Sc.newchanlabels{kk+1}=['ClassPr' num2str(kk)];
       Sc.newchantype{kk+1}='CLASSPR';     
    end;

    [ Dnew ] = osl_concat_spm_eeg_chans( Sc );    

    D.delete;
    D=Dnew;

    %%%%%%%%%%%%%%%%%%%
    %% generate source recon web report for this session
    report=osl_report_write(report);        
    source_recon_report=osl_report_add_sub_report(source_recon_report, report)    
 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% 

    source_recon_results.BF.data.D=spm_eeg_load(D);    
    source_recon_results.woi=woi;
    source_recon_results.samples2use=samples2use;
    source_recon_results.pca_order = pcadim;
    source_recon_results.normalisation = normalisation;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% save    
    
    source_recon_results.session_name=source_recon_sess.session_name;    
    source_recon_results.fname=[source_recon_results.session_name '_recon' ];
    disp(['Saving beamformer results: ' source_recon_results.fname]);
    
    osl_save_oat_results(oat,source_recon_results);

    results_fnames{sessi}=source_recon_results.fname;    

end

%%%%%%%%%%%%%%%%%%%
%% summary plots over sessions
source_recon_results.pca_order=nan(length(source_recon_results.pca_order),length(oat.source_recon.sessions_to_do),1);    
source_recon_results.normalisation=nan(length(source_recon_results.normalisation),length(oat.source_recon.sessions_to_do),1);    

oat.source_recon.results_fnames=results_fnames;

for sessi=1:length(oat.source_recon.sessions_to_do), sessnum=oat.source_recon.sessions_to_do(sessi);

    try,
        % load in opt results for this session:            
        res=osl_load_oat_results(oat, oat.source_recon.results_fnames{sessnum});

        mod_ind=find(strcmp(oat.source_recon.modalities,oat.first_level.report.modality_to_do));                
        
        for ff=1:length(res.pca_order),
            source_recon_results.pca_order(ff,sessi)=res.pca_order(ff);
        end;
        for ff=1:length(res.normalisation),
            source_recon_results.normalisation(ff,sessi)=res.normalisation(ff);
        end;

    catch ME,
        disp(['Could not get summary diagnostics for ' oat.source_recon.results_fnames{sessnum}]);
        ME.getReport
    end;
end;

for ff=1:length(res.pca_order),
    source_recon_report=osl_report_set_figs(source_recon_report,[source_recon_results.source_recon.modalities{ff} ' pca_order']);
    plot(oat.source_recon.sessions_to_do,source_recon_results.pca_order(ff,:),'*');xlabel('sess no.');ylabel([source_recon_results.source_recon.modalities{ff} ' PCA dim']); 
    source_recon_report=osl_report_print_figs(source_recon_report);
end;

for ff=1:length(res.normalisation),
    source_recon_report=osl_report_set_figs(source_recon_report,[source_recon_results.source_recon.modalities{ff} ' normalisation']);
    plot(oat.source_recon.sessions_to_do,source_recon_results.normalisation(ff,:),'*');xlabel('sess no.');ylabel([source_recon_results.source_recon.modalities{ff} ' normalisation']); 
    source_recon_report=osl_report_print_figs(source_recon_report);
end;

%%%%%%%%%%%%%%%%%%%
%% generate source recon web report
source_recon_report=osl_report_write(source_recon_report);        
source_recon_results.report=source_recon_report;

end
